Read data

source("scripts/myRLib.R")
funcset <- readLines("output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/func_subset.set")
funcset <- parseFuncSet(funcset)
funcset.names <- names(funcset)
xtyi.files <- strsplit("output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e+00.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p3.0000e-01.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e-01.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p3.0000e-02.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e-02.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p3.0000e-03.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e-03.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p3.0000e-04.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e-04.assoc.logistic", 
    ":")[[1]]
ps <- strsplit("1,0.3,0.1,0.03,0.01,0.003,0.001,0.0003,0.0001", 
    ",")[[1]]
inter.files <- strsplit("output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e+00.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p3.0000e-01.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e-01.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p3.0000e-02.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e-02.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p3.0000e-03.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e-03.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p3.0000e-04.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e-04.assoc.logistic", 
    ":")[[1]]
prs.files <- strsplit("../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e+00.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p3.0000e-01.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e-01.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p3.0000e-02.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e-02.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p3.0000e-03.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e-03.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p3.0000e-04.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e-04.txt", 
    ":")[[1]]

Epistasis

for (i in 1:length(xtyi.files)) {
    cat("##", funcset.names[i])
    cat("
")
    cat("
")
    snp.subset <- funcset[[funcset.names[i]]]
    
    # read I ~ SNP
    prs <- read.table(prs.files[i], header = T)
    prs$sid <- as.character(prs$sid)
    prs.sub <- prs[prs$sid %in% snp.subset, ]
    
    # read others
    xtyi <- read.table(xtyi.files[i], header = T)
    xtyi$SNP <- as.character(xtyi$SNP)
    yt.xtyi <- xtyi[(xtyi$SNP %in% snp.subset) & (xtyi$SNP %in% 
        prs.sub$sid), ]
    yt.xtyi.xt <- yt.xtyi[yt.xtyi$TEST == "ADD", ]
    yt.xtyi.yi <- yt.xtyi[yt.xtyi$TEST == "PRS", ]
    inter <- read.table(inter.files[i], header = T)
    inter$SNP <- as.character(inter$SNP)
    inter <- inter[inter$TEST == "ADDxPRS", ]
    yt.inter <- inter[(inter$SNP %in% snp.subset) & (inter$SNP %in% 
        prs.sub$sid), ]
    
    # order all rows
    yt.xtyi.xt <- yt.xtyi.xt[order(yt.xtyi.xt$SNP), ]
    yt.xtyi.yi <- yt.xtyi.yi[order(yt.xtyi.yi$SNP), ]
    yt.inter <- yt.inter[order(yt.inter$SNP), ]
    prs.sub <- prs.sub[order(prs.sub$sid), ]
    
    # check direction
    prs.sub$nt1 <- as.character(prs.sub$nt1)
    yt.xtyi.xt$A1 <- as.character(yt.xtyi.xt$A1)
    yt.inter$A1 <- as.character(yt.inter$A1)
    if (sum(prs.sub$nt1 != yt.xtyi.xt$A1) != 0 | sum(prs.sub$nt1 != 
        yt.inter$A1) != 0) {
        print("Direction does not match. Exit!")
        quit()
    }
    
    yt.xtyi.xt.mean.std <- getMeanStd(yt.xtyi.xt$SNP, yt.xtyi.xt$P, 
        yt.xtyi.xt$OR, "Y~Gi+I (Gi)")
    yt.xtyi.yi.mean.std <- getMeanStd(yt.xtyi.yi$SNP, yt.xtyi.yi$P, 
        yt.xtyi.yi$OR, "Y~Gi+I (I)")
    yt.inter.mean.std <- getMeanStd(yt.inter$SNP, yt.inter$P, 
        yt.inter$OR, "Y~Gi*I")
    prs.sub.mean.std <- data.frame(SNP = prs.sub$sid, Mean = prs.sub$ldpred_beta, 
        Std = 0, Type = "I~Gi")
    
    # create data frame for plot
    df <- rbind(yt.xtyi.xt.mean.std, yt.xtyi.yi.mean.std, prs.sub.mean.std, 
        yt.inter.mean.std)
    p <- ggplot(df) + geom_point(aes(x = SNP, y = Mean)) + geom_errorbar(aes(x = SNP, 
        ymin = Mean - 1.96 * Std, ymax = Mean + 1.96 * Std), 
        colour = "black", width = 0.1) + geom_hline(yintercept = 0, 
        color = "red") + facet_grid(. ~ Type, scales = "free_x") + 
        coord_flip() + theme(axis.text.x = element_text(angle = 90, 
        hjust = 1))
    p2 <- ggplot(yt.inter) + geom_bar(aes(x = SNP, y = -log10(P)), 
        stat = "identity") + coord_flip() + ggtitle("Interaction P-value")
    print(p)
    print(p2)
    cat("
")
    cat("
")
}

p1.0000e+00

p3.0000e-01

p1.0000e-01

p3.0000e-02

p1.0000e-02

p3.0000e-03

p1.0000e-03

p3.0000e-04

p1.0000e-04